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Abstract 

Explicit relations of matrices for two-dimensional finite element method with third-order trian- 
gular elements are given. They are more simple than relations presented in other works and could 
be easily implemented in new algorithms for both isotropic and anisotropic materials. Numeri- 
cal examples are given comparatively using second-order and third-order triangular elements for 
problems of wave propagation in rectangular waveguides which have analytic solutions. 
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I. INTRODUCTION 



The finite element method is a widely applicable numerical technique for obtaining ap- 
proximate solutions to boundary- value problems of mathematical physics [l, 2]. Any complex 
shape of the problem domain can be handed with ease by division into many subdomains, 
each subdomain being called a finite element. For two-dimensional problems one resorts usu- 
ally to triangular elements: the first-order triangular element, which requires three nodes, 
and the second-order triangular element, which requires six nodes. In order to achieve 
higher accuracy in the finite element solution, two approaches are commonly taken: one 
resorts to finer subdivision or smaller elements, and the other resorts to higher-order inter- 
polation functions or higher-order elements. Here we are interested in the two-dimensional 
finite element method with third-order triangular elements, each element requiring ten nodes 
which are numbered counterclockwise as shown in FigJU These high-order elements can be 
successfully employed for the characterization of wave propagation in shielded microstrip 
transmission lines or integrated circuits with slot lines, when the lines are infinitesimally 
thin and one needs to place a set of nodes above the line as well bellow the line, as if the 
line had a finite thickness Q . Generally, the elemental matrices for high-order elements are 
determined numerically, but a higher accuracy is assured with explicit expressions. Rela- 
tions for triangular elements, including the third-order ones, have been presented in [3]. In 
this paper the elemental matrices for third-order triangles are given in simple, explicit forms 
which can be easily implemented in different algorithms of waveguiding problems related to 
both isotropic and anisotropic materials. 




FIG. 1: Third-order triangular element 
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II. GENERAL RELATIONS 



Let us consider a two-dimensional boundary-value problem defined by the second-order 
differential equation 



dx\ x dx J dy\ v dy 



+ 



f, (x,y)eQ 



(1) 



where cj) is the unknown function; a x ,a y , and (3 are known parameters; and / is the source 
or excitation function. The ordinary two-dimensional Laplace equation, Poisson equation, 
and Helmholtz equation are special forms of (CQ). For simplicity we consider f = and the 
homogeneous Neumann boundary condition on the boundary enclosing the area Q. Within 
each element, <fi can be approximated as 

10 

<j> e (x,y) = '£Nfa,y)<l>° j , (2) 

where ^ are constant expansion coefficients and Nj(x, y) are the interpolation or expansion 
functions given by 
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In the above, the area coordinates L e - are given by 

[<ij I h r r ■ Cji/). 



2A e 



j = 1,2,3, 



(4) 
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in which cij , bj , and Cj are 



ai = x 2 y 3 - IJ2X3, 
a 2 = x 3 yx-y 3 xi, 
% = xxy 2 - yix 2 , 



h=y 2 - Va, 
h = Vs- Vx, 
h = Vi- V2, 



ci = x 3 - x 2 , 

C 2 =Xi- x 3 , 
C3 = X 2 — X\, 



where Xj and yj (j — 1, 2, 3) denote the coordinate values of the vertices and A e is the area 
of the eth element, 



-(bic 2 - b 2 ci) 



With the expansion of given in (j2J), we can proceed to formulate the elemental equations 



in pre 



using either the Ritz or Galerkin method [2|, |4J. As for example, in the Ritz method we 
formulate the problem in terms of a functional F((p) whose minimum corresponds to the 
differential equation of the boundary-value problem. The functional can be written as 



M 



(5) 



e=l 



where M denotes the total number of elements and F e is the subfunctional corresponding 
to the eth element. For differential equation (CQ), F e is defined as 



1 

2j.ni 



Oix\ 



fd(j) 



e x 2 



(6) 



V dx J ' y V dy 

where Q e is the domain of the eth element. Introducing expression (j2J) for e and differen- 



tiating with respect to 0f yields 



90| 
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j=x 



dN e dN e 8N e dN e 
ax ax oy ay J 



1,2,... 10. (7) 



In matrix form, this can be written as 



[IT] {0 e } 



(8) 



where 



dF e dF e dF e 



1 T 



with T denoting a transpose and the elements of the matrix [K e ] given by 



K e 



, a dNf aiv; 

n e " x dx dx 



+ a 



dNf ONI 
y dy dy 



+ (3N*m)dxdy, 



1,2,. ..10. 



(9) 
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III. ELEMENTAL MATRICES 



Assuming a x , a y , and /3 are constant in each element, we split [K e ] into three parts 

[K e ]=a x [Al]+a y [Al] + P[B% (10) 
where the elements of matrices [A%.], [Ay], and [B e ] are given by 

B% = JJ^N?N?dxdy, i, j = 1,2,... 10. (11) 

The integral calculus in (1111) can be performed analytically by using the following convenient 
integration formula for the area coordinates [sj 

J [jLty(L e 2 y(Ll) k dxdy = ( . + ^ ! +2) , 2A e , i,j, k = 0, 1, 2, 3, . . . (12) 

Explicit expressions for the matrix elements Bfj and A e xi ^ are given in Tables [I] and HT1 
respectively. Both [B e ] and [A*] are symmetric matrices, i.e., B? i = Bfj and A%^ — A%^. In 
Table HT1 we used notation 6^ = b^bj, with i, j = 1, 2, 3. The matrix elements • are obtained 
by changing 6^ to q (? = 1, 2, 3) in Table [TT1 

TABLE I: Matrix elements Bf- (i, j = 1, 2, . . . 10) are tabulated values multiplied by (3A e /2240) 
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TABLE II: Matrix elements A xi j = 1,2, ...10) are tabulated expressions multiplied by 

[81/(8A e )] 
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In the vector formulation of different waveguiding problems, the subfunctional F e of the 
eth element in ()6]) is more complicated and the matrix [K e ] in ffTUj) contains more elemental 
matrices. As for example, within a closed waveguide, the magnetic field satisfies the vector 
differential equation 



V x (—V x H) - = in fi 



(13) 



and the boundary condition n x (V x H) = on Ti, where Q denotes the cross section of 
the structure comprised by the electric wall Ti, ko is the wave number in vacuum, e r and 
fi r are the permittivity and permeability of the structure. With the z-dependence of H as 
H(x,y,z) = H(x, y)e^ UJt ~ kzZ \ where u is the circular frequency, the functional of (|T3|) can 

1 



be written as 
F(H 




+ jk z H y 



+ 



dy 

kfa r (\H x \ 2 + \H y \ 2 + \H z \ 2 )]d{l. 



jk z H x + 



dE 7 



dx 



+ 



dH v dH x 



dx 



dy 



(14) 



To render this as a real system, we introduce the transformation h z = —jH z , and with this, 
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pijl becomes 

F(H) 
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(15) 



The functional can readily be discretized in a standard manner, and the result is 
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(16) 



where the matrices are assembled from their corresponding elemental matrices, given by 

1 fd{N e }d{N e } T N 



JJtt e e r \ 



dy dy 

1 fd{N e }d{N e } T 
dx dx 



+ k 2 z {N e }{N e } T dtt 
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For simplicity we denote the matrices [C xy ] 
elements 



[Cy X ] , [D x ], and [Dy] having the matrix 
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dN ! dN L ^ 
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TABLE III: Matrix elements C% yi ^ (i,j = 1,2, ...10) are tabulated expressions multiplied by 
[81/(16A e )] 
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-P13-P31 

60 


-P13-P31 

60 


8 


-P21-9P31 


-P22 


18pi3"P23 


-P23-P32 


-P23-P32 
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Explicit expressions for the matrix elements C% yi j and are given in Tables [TTT] and HV] 
respectively. In Table HTT1 we used notation pij = biCj, with — 1,2,3. The matrix elements 
D^, are obtained by changing bi to q — 2, 3) in Table HVl 



S 



TABLE IV: Matrix elements = 1,2, ...10) are tabulated expressions multiplied by 

(27/140) 
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81 
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324 


19&3 
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72 
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561-863 
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2263-61 

72 


61 

12 


2 


196], 
324 
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81 
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562-861 
72 


2261-62 
72 


2263-62 
72 
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72 


-56 2 
72 


62 
12 
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196i 
324 


1962 
324 


I663 
81 


-56 3 
72 


-56 3 
72 


563-862 
72 


2262-63 
72 


2261-63 
72 


563-861 
72 


63 
12 


A 


23&i 
72 


-136 2 
72 


72 


-63 


6 2 -36i 
8 


261-63 
8 


61-63 
8 


62-263 
8 


63-62 
2 


263-361 
4 


5 


-13&1 
72 


236 2 
72 


56s 
72 


6i-36 2 
g 


-63 


63-61 
2 


61-263 
8 
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263-362 
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72 


-136 3 
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63-26! 
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-6l 


63-362 
8 


26 2 -6i 
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62-61 
8 


26i-36 2 
4 
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5&i 
72 
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236 3 
72 


63-61 
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263-61 
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-h 
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62-261 
8 


261-363 
4 
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-13&i 
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56 2 
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63-62 
8 


6i-26 2 
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62-61 
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-&2 


61-363 
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26 2 -3&3 
4 
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236i 
72 


56 2 
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-136 3 
72 


62—63 
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6 3 -26 2 
8 


61-62 
8 


26i-6 2 

8 


63-36i 
8 


-62 


26 2 -36i 
4 


10 


=&1 
12 


-&2 
12 


-63 
12 


361-263 
4 


362-263 
4 


36 2 -26i 
4 


36 3 -26i 
4 


863-262 
4 


36i-26 2 
4 






For a triangle with coordinates of vertices (1,1), (2,1), and (1,2), the area is A e = 
0.5 and the maximum absolute values on the columns of these elemental matrices vary in 
the intervals: Max|S?-| = 0.0057 - 0.1446, Max|A^| = - 4.05, Max|C^.| = - 2.025, 
Max|D^ | = — 0.2411. We can see that the elements of matrix [B e ] are very small. Thus, 
they could be more accurate than those obtained by performing the integrations numerically. 
Note that relations presented here for the matrix elements are more simple than those 
presented in |3|, |6|. 



IV. EXAMPLES 

In the following we consider some waveguiding problems which have analytic solutions. 
The most simple is the eigenvalue problem of a hollow square waveguide [7]. In Table 
IVl we give the results for the wavenumber fc (cm -1 ) obtained by using the vectorial fi- 
nite element method in terms of the magnetic field components [H x , H y , H z ], when the 
propagation constant on the z-axis direction, k z equals zero. An uniform, one-directional 
mesh is considered with second-order and third-order triangular elements in the domain 
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x, y G [—0.5cm 0.5cm]. For second-order triangles, the number of elements is 18 and the 
total number of nodes is 49. For third-order triangles, the corresponding numbers are 18 
and 100. Explicit relations of matrices for the second-order triangular elements were given 
in sl, Q]. As seen in Table |V] the results obtained using third-order elements have better 
accuracy than those obtained with the same number of second-order elements. 



TABLE V: Results for the hollow square waveguide 



Analytic 


Second Order 


Third Order 


Solution 


Triangles 


Triangles 


3.1416 


3.1438 


3.1416 


4.4429 


4.4523 


4.4431 


6.2832 


6.3451 


6.2852 



As another example we consider a rectangular waveguide partially filled with a dielectric 
[?| of relative permittivity e r = 6 as illustrated in Fig. |2j Results for /^(cm" 1 ) are given 
comparatively in Table [VTl at k z = and 1cm -1 , the rectangular domain of 1cm x 1cm being 
uniformly divided in 24 second-order elements with 63 nodes or 24 third-order elements with 
130 nodes. 




0.6 



1 



FIG. 2: Cross section of the inhomogeneous rectangular waveguide 



Finally we consider a rectangular waveguide with a 2 : 1 width to height ratio completely 
filled with a ferrite material characterized by a relative permittivity e r = 2 and a relative 
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TABLE VI: Results for the inhomogeneous rectangular waveguide 





Analytic 


Second Order 


Third Order 


k z 


Solution 


Triangles 


Triangles 





1.7666 


1.7681 


1.7666 




2.3053 


2.3076 


2.3053 




2.6779 


2.6875 


2.6779 


1 


1.8310 


1.8545 


1.8531 




2.3460 


2.3543 


2.3521 




2.7125 


2.7931 


2.7042 



TABLE VII: Results for the ferrite completely filled rectangular waveguide 



Analytic 


Second Order 


Third Order 


Solution 


Triangles 


Triangles 


0.6654 


0.6659 


0.6654 


1.3307 


1.3445 


1.3307 


1.9961 


1.9458 


1.9961 



permeability tensor /j, r given by [6|] 



Mr 



3 j'0.8 
1 
-j0.8 3 



The analytical solution for the wave number ko of the nth mode is 



k ,, 



3.0 
16.72 



7T 



r? - 



(18) 



(19) 



Results are given comparatively in Table rvTTI for ^(cm^ 1 ) of the first three modes, at k z = 0, 
when the rectangular domain of 2cm x 1cm is uniformly divided in 18 second-order elements 
with 49 nodes or in 18 third-order elements with 100 nodes. We can see that in the case of 
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third-order elements the results agree with analytical solutions very well. 
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